#ifndef MULTIROOTNEWTONRAPHSON_H_
#define MULTIROOTNEWTONRAPHSON_H_

// C++ includes
#include <vector>
#include <functional>
using namespace std;

// GSL includes
#include <gsl/gsl_multiroots.h>

class MultirootNewtonRaphson
{
public:
	typedef std::function<void(const vector<double>& x, vector<double>& f)>
	Function;

	typedef std::function<void(const vector<double>& x, vector<vector<double>>& J)>
	Jacobian;

	typedef std::function<void(const vector<double>& x, vector<double>& f, vector<vector<double>>& J)>
	FunctionJacobian;

public:
	MultirootNewtonRaphson();

	~MultirootNewtonRaphson();

	void SetNumEquations(unsigned numEquations);

	void SetFunction(const Function& function);

	void SetJacobian(const Jacobian& jacobian);

	void SetFunctionJacobian(const FunctionJacobian& functionJacobian);

	void SetTolerance(double tolerance);

	void SetMaxNumIterations(unsigned maxNumIterations);

	bool Iterate(vector<double>& solution);

	void Solve(vector<double>& solution);

private:
	static int GSL_Function(const gsl_vector* x, void* params, gsl_vector* f);

	static int GSL_Jacobian(const gsl_vector* x, void* params, gsl_matrix* J);

	static int GSL_FunctionJacobian(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J);

	struct Parameters
	{
		Function function;

		Jacobian jacobian;

		FunctionJacobian functionJacobian;
	};

private:
	gsl_multiroot_fdfsolver* solver;

	unsigned numEquations;

	Parameters parameters;

	double tolerance;

	unsigned maxNumIterations;
};

#endif /* MULTIROOTNEWTONRAPHSON_H_ */



//////////////////////////////////////////////////////////////////////////
// Template
//#ifndef MULTIROOTNEWTONRAPHSON_H_
//#define MULTIROOTNEWTONRAPHSON_H_
//
//// C++ includes
//#include <vector>
//#include <functional>
//using namespace std;
//
//// GSL includes
//#include <gsl/gsl_multiroots.h>
//
//template<class Parameter>
//class MultirootNewtonRaphson
//{
//public:
//	typedef std::function<void(const vector<double>& x, const Parameter& param, vector<double>& f)>
//	Function;
//
//	typedef std::function<void(const vector<double>& x, const Parameter& param, vector<vector<double>>& J)>
//	Jacobian;
//
//	typedef std::function<void(const vector<double>& x, const Parameter& param, vector<double>& f, vector<vector<double>>& J)>
//	FunctionJacobian;
//
//public:
//	MultirootNewtonRaphson();
//
//	void SetNumEquations(unsigned numEquations);
//
//	void SetFunction(const Function& function);
//
//	void SetJacobian(const Jacobian& jacobian);
//
//	void SetFunctionJacobian(const FunctionJacobian& functionJacobian);
//
//	void SetTolerance(double tolerance);
//
//	void SetMaxNumIterations(unsigned maxNumIterations);
//
//	void Iterate(vector<double>& solution); // todo
//
//	void Solve(vector<double>& solution, const Parameter& param);
//
//private:
//	static int GSL_Function(const gsl_vector* x, void* params, gsl_vector* f);
//
//	static int GSL_Jacobian(const gsl_vector* x, void* params, gsl_matrix* J);
//
//	static int GSL_FunctionJacobian(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J);
//
//private:
//	unsigned numEquations;
//
//	Function function;
//
//	Jacobian jacobian;
//
//	FunctionJacobian functionJacobian;
//
//	double tolerance;
//
//	unsigned maxNumIterations;
//};
//
//template<class Parameter>
//MultirootNewtonRaphson<Parameter>::MultirootNewtonRaphson() :
//tolerance(1.0E-6), maxNumIterations(1000)
//{}
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::SetNumEquations(unsigned numEquations)
//{
//	this->numEquations = numEquations;
//}
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::SetFunction(const Function& function)
//{
//	this->function = function;
//}
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::SetJacobian(const Jacobian& jacobian)
//{
//	this->jacobian = jacobian;
//}
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::SetFunctionJacobian(const FunctionJacobian& functionJacobian)
//{
//	this->functionJacobian = functionJacobian;
//}
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::SetTolerance(double tolerance)
//{
//	this->tolerance = tolerance;
//}
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::SetMaxNumIterations(unsigned maxNumIterations)
//{
//	this->maxNumIterations = maxNumIterations;
//}
//
//template<class Parameter>
//struct NewtonRaphsonParams
//{
//	typename MultirootNewtonRaphson<Parameter>::Function f;
//
//	typename MultirootNewtonRaphson<Parameter>::Jacobian J;
//};
//
//template<class Parameter>
//void MultirootNewtonRaphson<Parameter>::Solve(vector<double>& solution, const Parameter& param)
//{
//	const gsl_multiroot_fdfsolver_type* T = gsl_multiroot_fdfsolver_gnewton;
//
//	gsl_multiroot_fdfsolver* solver = gsl_multiroot_fdfsolver_alloc(T, numEquations);
//
//	gsl_multiroot_function_fdf func;
//
//	NewtonRaphsonParams<Parameter> params;
//
//	params.f = function;
//	params.J = jacobian;
//
//	func.n = numEquations;
//	func.f = &GSL_Function;
//	func.df = &GSL_Jacobian;
//	func.fdf = &GSL_FunctionJacobian;
//	func.params = &params;
//
//	gsl_vector* x = gsl_vector_alloc(numEquations);
//
//	for(unsigned i = 0; i < numEquations; ++i)
//		gsl_vector_set(x, i, solution[i]);
//
//	gsl_multiroot_fdfsolver_set(solver, &func, x);
//
//	int status = GSL_CONTINUE;
//
//	for(unsigned iter = 1; status == GSL_CONTINUE && iter <= maxNumIterations; ++iter)
//	{
//		status = gsl_multiroot_fdfsolver_iterate(solver);
//
//		status = gsl_multiroot_test_residual(solver->f, tolerance);
//	}
//
//	for(unsigned i = 0; i < numEquations; ++i)
//		solution[i] = gsl_vector_get(solver->x, i);
//}
//
//template<class Parameter>
//int MultirootNewtonRaphson<Parameter>::GSL_Function(const gsl_vector* x, void* params, gsl_vector* f)
//{
//	Function func = ((NewtonRaphsonParams<Parameter>*) params)->f;
//
//	vector<double> xv(x->data, x->data + x->size);
//
//	vector<double> fv(x->size);
//
//	func(xv, *(Parameter*)params, fv);
//
//	for(unsigned i = 0; i < x->size; ++i)
//		gsl_vector_set(f, i, fv[i]);
//
//	return GSL_SUCCESS;
//}
//
//template<class Parameter>
//int MultirootNewtonRaphson<Parameter>::GSL_Jacobian(const gsl_vector* x, void* params, gsl_matrix* J)
//{
//	Jacobian jac = ((NewtonRaphsonParams<Parameter>*) params)->J;
//
//	vector<double> xv(x->data, x->data + x->size);
//
//	vector<vector<double>> Jv(x->size, vector<double>(x->size));
//
//	jac(xv, *(Parameter*)params, Jv);
//
//	for(unsigned i = 0; i < x->size; ++i)
//		for(unsigned j = 0; j < x->size; ++j)
//			gsl_matrix_set(J, i, j, Jv[i][j]);
//
//	return GSL_SUCCESS;
//}
//
//template<class Parameter>
//int MultirootNewtonRaphson<Parameter>::GSL_FunctionJacobian(const gsl_vector* x, void* params, gsl_vector* f, gsl_matrix* J)
//{
//	GSL_Function(x, params, f);
//	GSL_Jacobian(x, params, J);
//
//	return GSL_SUCCESS;
//}
//
//#endif /* MULTIROOTNEWTONRAPHSON_H_ */
